Processed Splitting Algorithms for Rigid-Body Molecular Dynamics Simulations 
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A new approach for integration of motion in many-body systems of interacting polyatomic molecules is 
proposed. It is based on splitting time propagation of pseudo- variables in a modified phase space, while the real 
translational and orientational coordinates are decoded by processing transformations. This allows to overcome 
the barrier on the order of precision of the integration at a given number of force-torque evaluations per time 
step. Testing in dynamics of water versus previous methods shows that the obtained algorithms significantly 
improve the accuracy of the simulations without extra computational costs. 
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I. INTRODUCTION 

Systems of rigid bodies are widely used to model various 
phenomena on a broad range of length scales: from the mi- 
croscopic dynamics of molecules in gases and liquids (HQ], 
mesoscopic behavior of polymers and other complex collec- 
tions in chemical and biological physics [3, to macro- 
scopic movement of astrophysical objects in celestial me- 
chanics yila]. A lot of approaches, including the traditional 
Runge-Kutta and predictor-corrector schemes |l|] as well as 
more recent splitting techniques 10. HI S OH Ell EH]! have 
been devised over the years to integrate the rigid-body equa- 
tions of motion. 

Now it is well established that the most adequate integra- 
tion can be done by splitting the time propagator into analyt- 
ically solvable parts IT3I [l4i [l5ll . For Hamiltonian systems 
this provides the preservation of such essential properties as 
conservation of volume in phase space and time reversibility. 
As a result, the splitting algorithms exhibit remarkable stabil- 
ity and thus are ideal for long-duration molecular dynamics 
(MD) simulations. In addition, these algorithms can be sym- 
plectic, i.e. can exactly conserve the total energy associated 
with a nearby Hamiltonian. 

The splitting approach however has a limitation on the or- 
der K of precision at each given number n of force-torque 
evaluations per time step. Note that these evaluations present 
the most time-consuming part of the propagation. For this 
reason, the rigid-body motion in MD simulations is inte- 
grated mainly by the simplest ( K — 2) Verlet-type algorithms 
M |9l [lOl [Til 02] with n = 1. The optimized algorithms 
11511 at n = 2 can outperform Verlet schemes. But 
such an optimization does not rise the order of precision and 
for K = 2 only modest accuracy can be reached. Higher- 
order (K — 4) splitting schemes (note that K should be even 
to ensure time reversibility) can be derived beginning from 
n = 3 flU [l5tl . The grown computational costs at n — 3 and 
K — 4 can be compensated by the increased precision when 
adding gradient-like terms to the splitting propagator [ 15]. 

Meanwhile it has been found that the order K of precision 
can be risen by carrying out supplementary (so-called pro- 
cessed) decompositions apart from the basic (kernel) splitting 
(Trill . For K — 4, each minimal kernel and processor leads 
to one force and one force-gradient evaluations. This yields 



an effective number n = 2(1 + v), where v is the relative 
cost spent on the gradient evaluation with respect to that on 
the force calculation. This number can be decreased twice 
to n = 1 + v by constructing cheap approximate processors 
Taking into account that the evaluation of one 
force gradient is more expensive at least in a factor of v = 2 
than the calculation of one force B13lll4ll5ll gives that n > 3. 
However, the gradient evaluation may present a difficulty for 
systems with long-range (e.g. Coulomb) interactions, where 
the factor v can be too large IU5I1 because of the necessity to 
calculate cumbersome tail (Ewald-summated) contributions. 
Note also that the processed algorithms of Refs. lfl7l[l8ll were 
obtained exclusively for pure translational motion and they are 
not suitable for rigid-body dynamics. The processing meth- 
ods introduced in Refs. Ilia, 1 1 911 for solving ordinary differ- 
ential equations are more general but need an adaptation to 
be exploited in the case of rotational motion. In particular, 
contrary to free translational dynamics, the propagator of free 
rotational motion cannot be handled at once and requires ad- 
ditional splitting into analytically integrable parts [15] or in- 
volving special functions 12011 . 

Up to now, no processing schemes were designed and ap- 
plied to MD simulations of interacting rigid bodies. The 
rotational motion is much more complicated than transla- 
tional displacements and thus demands a separate investiga- 
tion. Moreover, a fundamental theoretical problem on the 
possibility to overcome the barrier n — 3 for the fourth-order 
integration still remains open. This overcoming is important 
from the practical point of view as well, because smaller val- 
ues of n could noticeably speed up the calculations in view of 
the restricted capabilities of even supercomputers. 

In the proposing paper we develop the processing formal- 
ism in the explicit presence of translational and orientational 
degrees of freedom. We show that using a proper transfor- 
mation of phase coordinates allows to lower the fourth-order 
barrier to the value n = 2 with no gradient evaluations. It is 
proven also that in a specific case of quasi-fourth-order inte- 
gration the number of force-gradient evaluations per step can 
be reduced to n = 1 at all. 

The paper is organized as follows. The new processed algo- 
rithms are consistently derived in Sec. II. Their applications 
to rigid-body MD simulations and comparison with integra- 
tors known previously are presented in Sec. III. Concluding 
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remarks are highlighted in Sec. IV. 



n. THEORY 

Let us consider a classical system of TV interacting rigid 
polyatomic molecules. The dynamical state of such a sys- 
tem in the laboratory frame is determined by the position 
Yi of the center of mass m of the zth molecule, its atti- 
tude matrix Sj as well as the translational p; and angular 
q ; momenta. The equations of motion can be written in 
the following compact form dp/dt = Lp(t). Here p — 
{ri,Pi,S 1 ,q 1 ;...;r A r,p A r,S A r,q A r} ee {r,p,S,q} is the 
set of phase variables, 
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f i (r,S)~ + g i (r,S)~ 
opi dqi 



(1) 



denotes the Liouville operator, fi and gj are the force and 
torque, respectively, acting on the molecule due to atomic in- 
teractions, 

tt z -fly 

w(n) = ( -n z o n x 
n Y -six o 

is the skewsymmetric matrix related to the principal compo- 
nents (tlx, fiyj &z) of the angular velocity Jl = J^Sq with 
J = diag( Jx, Jy, Jz) being the matrix of moments of iner- 
tia. If an initial configuration p(0) is specified, the unique 
solution to the equations of motion can formally be cast for 
any time t as p(t) — [exp(Lh)] k p(0), where h = t/k is the 
size of the time step and k denotes the total number of steps. 

In the standard splitting approach lfl3l \\A Il5ll . the Li- 
ouville operator L = A + B is decomposed into its ki- 
netic A = m- 1 p-d/dr + W(Ti)S-<9/<9S and potential B = 
f(r, S)-<9/<9p + g(r, S)-9/9q parts (we will omit the sub- 
script i for the sake of simplicity). Then the one-step time 
propagator e Lh can be factorized as e ( A + B ) h +°( h ) = 
n^ii e Bb » h e Aa » h ee $ K (h), where n > 1 and {a^ b^} are 
chosen in such a way to provide the highest possible order 
K of precision, and 0{h K+1 ) denotes the local error. For 
instance, the second-order (K = 2) Verlet algorithm is ob- 
tained at n = 1 by e (A+B)h+a(h 3 ) = ^^Ah^ = 
Note that the decomposition constants a M and 6 M should en- 
ter symmetrically into the factorization to ensure its time re- 
versibility. This reduces the total number of independent 
constants from 2(n + 1) to n + 1. In turn the symme- 
try provides automatic cancellation of all even-order terms 
in 0(h K+1 ), leading to evenness of K. For even orders 
K > 2, the local error function has the form 0(h K+1 ) = 
ci[A [A,B]]h 3 + c 2 [B, [A,B]]h 3 + 0(h K+3 ), where [ , ] 
designates the commutator operation and the coefficients c\ 
and C2 depend on {a^,^}. At K = 2, the two conditions 
^2 a p = bp = 1 should be satisfied to exclude the 
zeroth-order term from 0(h K+1 ). In order to increase the 



precision to K = 4 we should satisfy the two additional con- 
ditions ci({a M , b^}) = C2({a M , b^}) — 0. This can be pro- 
vided by increasing the number n + 1 of independent con- 
stants at least to the number of the order conditions, i.e, to 4. 
We see thus that fourth-order (K — 4) schemes can be con- 
structed only beginning from n = 3 and this number cannot 
be lowered within the standard splitting method. At n = 3, 
the fourth-order (K = 4) factorization can be presented as 
the concatenation $ 4 (/i) = $ 2 (x/i)3>2(l - 2x)h)<& 2 {xh) + 
0(h 5 ) of three Verlet signatures, where x — 1/(2 — v2). 

For arbitrary times t, the solution to the equations of mo- 
tion can be evaluated by consecutively applying k times the 
one-step splitting propagation ^x(h). This yields p(t) = 
[$K (h)] k p(0) + 0{h K ), where 0(h K ) ~ kO(h K+1 ) is the 
global error due to the accumulation of the local one after 
k = t/h 1 steps. The action of the exponential operators 
e Ar and e Br on a phase space point p is given analytically by 



e T {r,p,S,q} = {r + m 1 pr, p, S(q, r)S, q} , 
e BT {r, p, S, q} = {r, p + f (r, S)r, S, q + g(r, S)r}, 



(2) 



were the shift of r corresponds to free translational motion (at 
constant p), while the changes in p and q relate to motion in 
instantaneous force-torque fields 1151 . The matrix S(q, r) ex- 
actly propagates S over time r according to the free rotational 
dynamics (q remains constant) dS/dt = W(J _1 Sq)S. Ex- 
pressions for S(q, t) in terms of efficient routines for elliptic 
and theta functions are reported in Ref. [20]. Alternatively, 
H(q, r) can be replaced by its second- or fourth-order coun- 
terparts H 2 (r) = * x (§)*y(§)* z (r)*y(§)* x (§) and 
3 4 (t) = S 2 (xr)S 2 ((l - 2x)r)H 2 (xr), where * c (r) = 
cxp[W(f^)r] ee 0(fif,r) is the matrix representing rota- 
tion on angle Q^t around axis £ at constant component f\ of 
fi = J _1 Sq (see Eq. (19) of Ref. (H for 0(fi c ,r)). Note 
that each force-torque recalculation in e Br requires oc N 2 op- 
erations that is the most time-taking part of the splitting prop- 
agation, while the costs for handling e Ar are negligible (pro- 
portional to N). The total number of force-torque recalcula- 
tions per step in $a' is equal to n. 

The commutators [A, [A,B]] and [B, [A,B]] which appear 
in the local error function 0(h K+1 ) can be calculated explic- 
itly using the expressions for operators A and B. Then, in the 
case of the Verlet algorithm (K = 2) we find c\ = 1/12 = 
2c 2 and 0{h 3 ) = -{2rrr 1 i ■ d/dr - f • d/dp)h 3 /12 + 
0(h 5 ), where at the moment the orientational degrees of free- 
dom were frozen to simplify notation. Transferring now the 
corresponding parts of 0(h 3 ) from e Lh +°( h ) to the right 
under the exponentials e Ah and e B i one obtains e Lh = 
e B ie Ah e B i + 0{h 5 ), where A = A + m^f ■ d/drh 2 /6 
and B = B — f • d/dph 2 /12 are the modified counter- 
parts of A and B. Thus, the order of the Verlet signature 
can increase from K = 2 to K = 4 when the decomposi- 
tion is performed for the nearby Liouvillian C = A + B = 



l f • d/drh 2 /6 - f • d/dph 2 /12), where the equal- 



ities f = d{ jdt = Lf and f = Lf for the time derivatives of f 
have been applied. Note however that the nearby exponentials 
e Ar and e Br cannot be handled analytically in p-space (unlike 
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e and e , see Eq. (2)), because of the existence of compli- 
cated functions f = f (p) and f = f (p) which contrary to the 
force field f (r, S) depend not only on the positions (r, S) but 
on the momenta (p, q) as well. 

The main idea of our approach consists in finding such 
a processing transformation p = Xp from the phase space 
point p to a new set p of variables to make the action of 
the nearby exponentials analytically calculable. Taking into 
account the explicit structure for the nearby Liouvillian C, 
the general form of the desired transformation reads 1 = 
(r + am-Hh 2 )d/dv + (p + (3ih 2 )d/dp + 0(h 4 ) = 1 a ,p, 
where a and j3 are some coefficients which will be defined 
below. It can be verified readily that in the new variables, the 
equations of motion become dp/dt = Lp, where L = A + B 
is the corresponding Liouville operator with A = m _1 [p + 
(a - (3)f(r)h 2 ]-d/dr and B = [f(r - am^i^h 2 ) + 
/3i (?) h 2 ] -d/dp. Then for the nearby counterparts of A and B 
one finds A = A + m^f (?) • d/drh 2 /6 and B = B - f (?) • 
d/dph 2 /12. We see that the terms with f and f can be killed 
in A and B by putting (a-/3) = -l/6and/3 = 1/12, i.e. a = 
— 1/12. The orientational degrees of freedom can be included 
in a similar manner leading to the total processing transfor- 
mation 1 a . p = (r + am'Hh 2 )d/dr + (p + (iih 2 )d/dp + 
©(j^Sglr, S), ah 2 )Sd/dS+{q+[3gh 2 )d/d(i+0(h 4 ) and 
the nearby operators A = TO~ 1 p-9/9r + W(J~ 1 Sq)S-<9/<9S 
and B = f (? 7 , S 7 )-<9/<9p + g(? 7 , S 7 )-<9/<9q = B^ 
at a = -1/12, = 1/12, and 7 = -a = 
1/12. Here ©(j-^g^, S), ah 2 ) = cxp[W(Q)ah 2 } 
is the matrix representing three-dimensional rotation (i.e. 
0(Q,r) = Icos(Qr) + [1 - cos(Qr)] [W(Q)W(Q)/Q 2 + 
I] + sin(Qr)W(Q)/(5 with I being the unit matrix) around 
vector Q = J _1 Sg on angle Qah 2 , and {r 7 ,S 7 } — 
X 7i o{?, S} are the auxiliary position and attitude matrix. 

From the aforesaid, we have for the one-step propagation 
in p-space that p(t + h) = c Lh p(t) = e B -'ie Ah e B '>i p(t) + 
0(h 5 ). In p-space the solution can be reproduced by applying 
the inverse transformation 2~^g as pit + h) = c Lh p(t) — 
l^pP^ + h). This leads to the resulting propagation of p in 
the form 

e Lh = 1-^e B ^e Ah e B ^1 at0 + 0(h 5 ), (3) 

where a = -1/12, (3 = 1/12, and 7 = 1/12. The operator 
1 a , /3 transforms a phase space point p to the set p = 1 a ,pP = 
{?, p, S, q} of time-step dependent pseudo-variables, where 

r = r + am- 1 f(r,S)/i 2 , p = p + f3f (p)h 2 , 

(4) 

S = ©(j^Sgtr, S), ah 2 )S, q = q + /3g(p)/i 2 . 

The action of the exponential operators c At and c B ~< T can be 
given analytically as 

c At p = {r + m _1 pr, p, H(q, r)S, q}, 

e B - T p = {?, p + f(r 7 , S 7 )r, S, q + g(? 7 , S 7 )r}. 

Expressions (5) are similar to Eq. (2), since besides the formal 
replacement of p by p the only difference between (A, B) 



and (A,B) lies in the modification of the force f(r 7 ,S 7 ) 
and torque g(? 7 ,S 7 ). Apart from the calculation of their 
basic values f(r, S) and g(r, S), the modification requires 
(for 7 /: 0) one extra force-torque evaluation at the auxil- 
iary positional r 7 = r + 7m~ 1 f (?, S)h 2 and orientational 
S 7 = ©(j _1 Sg(r, S), "fh 2 )S coordinates. This increases 

the number of force-torque calculations in c B '< T from n = 1 
(at 7 = 0) to n = 2 (at 7 7^ 0), but the order of precision 
of the processed splitting propagation grows from K = 2 (at 
a — (3 = 7 = when it reduces to the genuine Verlet signa- 
ture) to K = 4 (at -a = (3 = 7 = 1/12). 

Because of T~^T Qj /3 = 1, the solution to the equa- 
tions of motion can now be cast for any t as p(t) = 

e K]p\ e ^ e Ah e 6 ^] kc Z a ^p{Q) + 0(h 4 ). Then the process- 
ing transformation X a>/ g can be performed only once on the 
very beginning, while the inverse transformation only 
once at the end of the considered time interval [0,i\. In view 
of this, the step by step integration can be interpreted as the 
time propagation of pseudo-variables p by the kernel splitting 
e B i 2 Q Ah Q B ~i 2 in the transformed phase space. The real phase 
coordinates p are not involved explicitly into the consecutive 
updating process. They can be reproduced from p whenever it 
is necessary (for example, when the measurement is desired) 
using the inverse transformation p = T~^p. This transfor- 
mation reads (cf. to Eq. (4)) 

r = f-am- 1 f(f,S)/i 2 , p = p - f3f(p)h 2 , 

S = 0(-J- 1 Sg(r,S), a / i 2 )S, q = q-/3g(p> 2 , (6) 

where the higher-order terms 0(h 4 ) have been neglected 
since they are not accumulated in p(t). 

The next crucial point concerns the evaluation of time 
derivatives f(p) and g(p) which arise in Eq. (6). It is ob- 
vious that their direct evaluation should be obviated since 
this results in complicated gradient terms. Fortunately, the 
derivatives can be evaluated at a given t in a quite efficient 
way by the symmetric interpolation {f , g}(p) = [{f , g}(t + 
h) - {f , g}(i - h)]/(2h) + 0{h 2 ), where {f , g}(i ± h) = 
{f, g}(r(t±h), S(t±h)). Such an interpolation is indeed re- 
alizable because the pseudo-variables p(t±h) are determined 
step by step in the course of the kernel propagation indepen- 
dently of p(t). Then the real variables p{t) can be reproduced 
from p(t) with a one-step retardation, when the pseudo-phase 
coordinates were already propagated to p(t + h). This avoids 
the calculation of extra forces and torques during the inter- 
polation and involves only those which already were eval- 
uated within the kernel propagation. The time derivatives 
f(p) and g(p) in Eq. (4) can be evaluated as {f, g}(p) = 
[{f,g}(f ) - {f,g}(-f )]//* + 0(h 2 ), where {f,g}(±f ) = 
{f,g}(r(±f),S(±f)) withr(±|) = r(0)±m^ 1 p(0)| and 
S(±f) = efiJ-^OMO), |)S(0). This involves two ex- 
tra forces and torques at ±h/2 but exclusively on the first step 
of the integration when starting from an initial configuration 
p(0) and performing the direct transformation 1 a ,p- 

We see therefore that the processed splitting (PS) algorithm 
derived is truly of the fourth order and requires only n = 2 
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force-torque evaluations per time step. This overcomes the 
barrier n — 3 inherent in standard schemes. Moreover, the al- 
gorithm is time reversible [because the exponential operators 
enter symmetrically into the propagator (Eq. (3))] and phase- 
area preserving [since simple shifts and rotations (Eq. (5)) do 
not change the volume]. In addition, the algorithm is explicit 
(no iterations) and exactly conserves the rigid molecular struc- 
ture (because H and are rotational matrices). The kernel 
splitting can also be made symplectic, because it is based on 
the Verlet-like signature which at 7 = conserves a nearby 
Hamiltonian II 15ll20l 12111 ). For a finite 7^0, the potential op- 
erator can be represented by B 1 = Bo + 7[£?o, [A, Bo]]h 2 /2 + 
0(h 4 ), where [B ,[A,B ]] = (B 6 - B )/e + 0(eh 4 ) 
with e <C 1. Then the modified force and torque in B 1 
can be evaluated as f(f 7 ,S 7 ) = f(r, S) + 7Af(r, S) and 
g(f 7 ,S 7 ) = g(f , S) + 7Ag(r, S), where the secondary 
fields are Af(f , S) = [f(r e , S e ) - f(r, S)]/e + 0{eh 4 ) and 
Ag(r, S) - [g(r E , S e ) - g(r, S)]/e + 0(eh 4 ). The param- 
eter e is typically taken to be of order 10~ 4 for double preci- 
sion arithmetic to minimize the effect of 0(eh 4 ) -terms while 
avoiding round-off truncations. The processing transforma- 
tions (Eqs. (4) and (6)) need not be necessarily symplectic, 
since their effects are not propagated (2T^ g^a,0 = !)■ 

That is very surprising, within the PS method the num- 
ber n of force-torque recalculations per time step can be 
reduced to n = 1 at all when a quasi-fourth order is re- 
quested. Note that the true fourth order means that the de- 
viations of the generated trajectories p(t) from their exact 
counterparts are equal to 0{h 4 ) ~ Ch 4 at t » h. In MD 
simulations, this strong requirement may not be so needed, 
because according to the Lyapunov theorem [3] the coeffi- 
cient C ~ e At grows (A > 0) exponentially with increas- 
ing t. Then the concept of the quasi-fourth order can be 
more useful. It implies that the deviations apply not to in- 
dividual variables of each particle but rather to a collective 
function for which C is independent of t. In microcanonical 
simulations such a function should be the total energy E = 
\ EZM 2 /rn + ClM) + \ £g^ 6 ^(|r 4 - Tj |, S,, Sj) 
of the system, where ip ab denotes the intermolecular atom- 
atom potentials, and M is the number of atoms per molecule. 
Cumbersome analysis shows that E can be conserved with 
the fourth-order accuracy at ji = 1 by tuning the parame- 
ters of the method to a = -1/24, [3 = 1/12, 7 = 0, and 
rj = 1 /48 (then p(t) and other quantities will not be necessar- 
ily reproduced up to the fourth order). Here we should add a 
new 77-term when transforming (Eq. (6)) angular momentum 
asq = q-/3g(p)/i 2 +7 ? S+[JW(J- 1 Sq)J- 1 -W(J- 1 Sq)- 
W(Sq)J _1 ]Sg(r , S)h 2 , where S + is the transposed matrix 
(and correspondingly modify Eq. (4)). This adding presents 
no difficulty since g(f , S) was already calculated during the 
kernel splitting. For systems without periodic boundary con- 
ditions, e.g. in celestial mechanics, the total angular momen- 
tum is often also conserved. It will be kept with the second- 
order accuracy by the quasi-fourth integrator (n = 1). This 
is in contrast to the genuine fourth-order algorithm (n = 2) 
which produces all quantities to within the 0(h 4 ) precision. 
Therefore, the former integrator may be less universally appli- 



cable than the latter one. The PS algorithms will be referred 
to as PS1 (n = 1) and PS2 (n = 2), respectively. 

Further improvements are possible by splitting the atom- 
atom potentials into short- and long-range parts. Then a 
multiple-time stepping (MTS) technique M22I1 can be em- 
ployed, where the expensive long-range (weak) forces are 
sampled less frequently using larger time steps, while the 
short-range (strong) interactions are integrated more accu- 
rately inside the kernel propagator using smaller steps. The 
MTS implementation within the PS method goes beyond the 
scope of this paper and will be considered elsewhere. 

III. NUMERICAL RESULTS 

We first present the proposed PS method (see Sec. II) 
in algorithmic form to simplify its implementation in a 
numerical code. Thus, starting from an initial configu- 
ration p(0) = {r(0),p(0),S(0),q(0)} at t = and 
calculating the three forces f(0) and f(±/i/2) as well 
as the three torques g(0) and g(±/i/2) at the positions 
{r(0),S(0)} and {r(±/i/2), S(/i/2)}, respectively, where 
r{±h/2) = r(0) ± m _1 p(0)/i/2 and S(±h/2) = 0( ± 
J _1 S(0)q(0), /i/2)S(0), we make the direct processing 

transformation (Eq. (4)) to p(0) = {r(0), p(0), S(0), q(0)} 
as 

f(0)=r(0)+am" 1 f(0)/i 2 , 

S(0) = ©(J- 1 S(0)g(0),^ 2 )S(0), 
p(0)=p(0)+/3(f(ty2)-f(ty2))fc, 

q(0)=q(0)+/3(g(V2)-g(V2))/i- 

Having p(0), we calculate the two initial forces f(0) and 
f e (0) as well as the two initial torques g(0) and g E (0) 
at the positions {f(0),S(0)} and {f e (0), S £ (0)}, respec- 
tively, where r s (0) = r(0) + £m -1 f(0)/i 2 and S £ (0) = 
0(J- 1 S(O)g(O), eh 2 )S(Q). Note that the direct transforma- 
tion (Eq. (7)) as well as the evaluation of the initial forces and 
torques should be carried out only once at the very beginning 
(t = 0) of the integration. 

Now we perform the single-step propagations of p from 
time t to t + h according to the kernel splitting (Eq. (3)) as 
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where P t +^ and qt+^ are the intermediate values, and the 
two new forces f (t + h) and f e (t + h) as well as the two new 
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torques g(i + h) and g E (t + ft) should be calculated at the 
new positions {r (t + ft) , S (t + ft) } and {r E (t + h), S £ (t + h)}, 
respectively, with r e (t + ft) = f (t + h) + em^ 1 i(t + h)h 2 and 
S e {t+h) = 0(J- l S(t+/i)g(t+/i),£/i 2 )S(t+/i) before the 
evaluation of p(t+/i) andq(i + ft). Saving the forces f(t+ ft) 
and f e (f + ft,) as well as the torques g(t + ft) and g E (t + ft), 
we repeat Eq. (8) (with formal replacing t by t + ft in it) to 
propagate p from time t + ft to t + 2ft. In such a way, step by 
step we can recycle Eq. (8) arbitrarily number k > 1 of times 
and obtain the value of p(t) for any t = kh. Each recycle will 
require the recalculation of only two (n = 2) new forces and 
torques. 

When at least two recycles of Eq. (8) are done already, we 
will have the three consecutive values pit — ft), pit), and 
pit + ft) for some t — kh. The forces f(i) and f(i ± ft) 
as well as the torques g(i) and g(t ± ft) will also be already 
known because of the kernel propagations. Then we can make 
the inverse processing transformation (Eq. (6)) of pit) to the 
genuine value pit) at a current t according to 

r(i) = f{t) - am-^i^h 2 , 
S(t) = 0(-J- 1 S(t)g(t), a ft 2 )S(i), 
pit) = pit) - 0(f (t + ft) - f (t - ft)) ft/2, 
q(i) = q(i) - /3(g(i + ft) - g(f - ft)) ft/2, 

and calculate at this point all necessary observable quantities 
(such as the total energy, etc.). This completes the PS2 algo- 
rithm in = 2), where a = -1/12, j3 = 1/12, and 7 = 1/12. 
ThePSl integrator^ = 1) follows at a = -1/24, (3= 1/12, 
7 = 0, and 77 = 1 /48 (here the evaluation of the modified 
force f e and torque g e should be omitted in Eq. (8) since 
7 = 0, while the inclusion of the ?7-term in Eqs. (7) and (9) is 
trivial). 

For testing of the algorithms we applied the TIP4P model 
(M = 4) of water Hf with N — 512 molecules. The MD 
simulations were carried in the microcanonical (NVE) en- 
semble at a density of N/V — 1 g/cm 3 and a temperature of 
292 K. The Ewald summation 124ll was exploited to handle 
long-range Coulombic atom interactions. The accuracy of the 
simulations was measured by calculating the ratio 1Z of the 
fluctuations of the total energy E to the fluctuations of its po- 
tential part [15]. The computational costs T were estimated 
in terms of the number of force-torque evaluations in a given 
time interval, taken to be A = 1 ps, so that T = nA/h. The 
equations of motion were solved at several sizes of the time 
step ranging from ft = 0.5 fs to 5 fs. In total k = t/h = 10 5 
steps were used for each algorithm and each step size. 

The costs T versus precisions 7Z of the integration obtained 
within the two proposed PS algorithms (K — 4) at the end 
of the simulations are plotted in Fig. 1 by the curves marked 
as PS1 in = 1) and PS2 (n = 2), respectively. The results 
corresponding to the Verlet-type (VT) algorithm (K = 2 and 
n = 1), its optimized (VO) version (AT = 2 and n = 2), the 
Forest-Ruth (FR) scheme iK = 4 and n = 3), as well as the 
gradient-like (GL) algorithm iK = 4 and n — 3) (these inte- 
grators are described in Ref. 11151 l2lll ) were also included for 
the purpose of comparison. It has been established that other 
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FIG. 1 : The cost versus relative error for different algorithms in MD 
simulations of water. The circles correspond to the time steps (left to 
right) h = 1, 2, 3, 4, and 5 fs. The dashed lines represent the most 
characteristic levels. 

known rigid-body integrators fiiMEIlIlllII iK = 2 and 
n = 1) behave similarly to the VT algorithm. Higher-order 
schemes lfl5l l26ll with K > 4 and n > 4 are less efficient in 
MD simulations because of the large numbers of costly force- 
torque recalculations. The processed fourth-order algorithm 
by Blanes and Casas (BC) et al. d [Hi with K = 4 and 
n = 1 + v = 3 (where the kernel and processor are defined 
according to Eqs. (20) and (21) of Ref. [16]) was adapted to 
rigid-body motion and considered too. 

As can be seen from Fig. 1, with decreasing T (rising ft) 
each curve terminates at some point where the simulations 
begin to exhibit a drift in 1Z. This happens around ft ~ 5 
fs (larger ft can be used within the MTS). At the minimally 
possible costs T ~ 200, the VT integrator can provide only 
a crude energy conservation 7Z ~ 7%. This level of errors is 
too large and generally unacceptable in MD simulations. It 
should be reduced at least to 1Z ~ 1%, arguably the upper 
limit of allowable error for which the dynamics can be sim- 
ulated adequately. The proposed PS 1 algorithm just satisfies 
this criteria even at T ~ 200. On the other hand, the level 
1Z ~ 1% can be achieved by the VT integrator by increas- 
ing the load to T ~ 550, i.e. in a factor of 2.75. Thus the 
PS 1 algorithm may spend considerably smaller CPU time at 
a given precision. The PS2 algorithm is also superior to the 
VT scheme. For more accurate (1Z < 1%) simulations, the 
relative efficient of the PS algorithms (AT = 4) with respect 
to the VT scheme (AT = 2) rises further (because 1Z ~ ft A ) 
and reaches a factor of 5 at 1Z ~ 0.1%. At the same time, for 
T ~ 550 the PS2 and PS 1 algorithms are able to lower the 
numerical errors from the value 1Z ~ 1% inherent in the VT 
integrator to the levels 1Z ~ 0.2% and 0.02%, respectively, i.e. 
up in 50 times! The VO integrator is clearly inferior to the PS 
algorithms, although it can be better than the VT signature. 
The BC scheme is superior to the VO integrator but worse 
than the PS algorithms. The FR scheme leads to the worst ef- 
ficiency. The GL algorithm can be used only at T > 750, i.e. 
when a very high accuracy (1Z < 0.02%) is required. Then it 
appears to be more efficient than the PS2 integrator. However, 
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FIG. 2: The fluctuations (a) and deviations (b) of the total energy 
versus the length of the MD simulations carried out at h — 4 fs using 
different algorithms. 

the PS 1 algorithm is the best in the whole T-region. 

Samples of the relative fluctuations lZ(t) and normalized 
deviations SE(t) = (E(t) -E(0))/E(Q) of the instantaneous 
total energy E(t) are shown in subsets (a) and (b) of Fig. 2, re- 
spectively, versus the length t/h of the simulations performed 
at a typical step h = 4 fs using different integrators. We can 
observe in Fig. 2(a) that the functions lZ(t) are flat with no 
drift on the entire time domain. The PS algorithms thus apart 
from their high efficiency, exhibit also excellent stability prop- 
erties. As is illustrated in Fig. 2(b) for the PS1 method, the 
total energy E{t) continues to keep near its initial value E(0) 
even after an extremely long period of time with k = 10 6 



steps. The magnitude of the deviations 8E(t) is quite small 
and does not exceed a level of 0.01%, making the energy con- 
servation almost exact. 



IV. CONCLUSION 

In this paper we have proposed a novel method for the 
integration of motion in rigid-body MD simulations. It 
combines standard splitting techniques with special phase- 
space processing transformations. Comparison with the well- 
recognized previous schemes has demonstrated that the new 
method allows to significantly improve the efficiency of the 
integration with no extra computational costs. The algorithms 
obtained are easy in implementation and can readily be incor- 
porated into existing MD codes. They can also be applied to 
hybrid Monte-Carlo, MD simulations of simple fluids and to 
other fields mentioned in the introduction as well as be ex- 
tended to more complicated systems with flexible molecules. 
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